*Test rank invariance/similarity at individual quantiles using the JTPA training data - Table 3, Panel A: earnings, I: Individual Quantile Test
*Please make sure `ivqte' and `moremata' are installed
cd "C:\Users\XXX"
*Replace `Gender' by male or female
log using ".\Dropbox\Rank Invariance\Empirical application\JTPA\\Results\Table3_IA_Gender.smcl", replace
clear all
set more off
set matsize 5000
version 12

	capture program drop test_indqt
	program define test_indqt, eclass 
	tempname diff qte qt_c qt_t
	tempvar qtindicator
	matrix `diff'=J(1,dfV,.)
	tempvar y2
	qui gen `y2'=y
	qui replace `y2'=0 if t==1
	qui ivqte y (t=z), quantiles($qt) dummy($cov)
	matrix `qte'=e(b)
	qui ivqte `y2' (t=z), quantiles($qt) dummy($cov)   // use y2 to trick ivqte estimate quantile under control
	matrix `qt_c'=-e(b) // notice the negative sign
	matrix `qt_t'=-e(b)+`qte'
	local qind=0
	gen `qtindicator'=0
    foreach indqt of global qt {
	local qind=`qind'+1
	qui replace `qtindicator'=(y<=`qt_c'[1,`qind']) if t==0
	qui replace `qtindicator'=(y<=`qt_t'[1,`qind']) if t==1
	local lind=0
	foreach l of global level {
		local lind=`lind'+1
		qui sum `qtindicator' if z==1 & x==`l'
		scalar temp=r(mean)
	    qui sum `qtindicator' if z==0 & x==`l'
		matrix `diff'[1,numqt*(`lind'-1)+`qind']= r(mean)-temp
		}
	}
	ereturn post `diff'
	end 

use ".\Dropbox\Rank Invariance\Empirical application\JTPA\JTPA_ContinousAge.dta", clear
global qt "0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85"
gen y=earn
gen z=assignmt
gen t=training
qui sum y
local obs=r(N)
replace hsorged=round(hsorged)
replace married=round(married)
replace wkless13=round(wkless13)
*Table 3: Female
local male=0
**Table 3: Male
*local male=1

*cov =0/1 controls for whether using covariates in the first-stage QTE estimation.
local cov=1
local mincell=5
local bsreps=2000

*age is in 6 categories, which are represented by 5 age dummies
local cov1 "hsorged  black   hispanic  married  wkless13  age2225 age2629   age3035  age3644   age4554 "
local cov0 "hsorged   black    hispanic   married   wkless13 afdc age2225    age2629    age3035 age3644   age4554 "
qui tostring `cov1', replace
qui tostring `cov0', replace
*sex==1 for male, sex==0 for female	
gen x=hsorged+black+hispanic+married+wkless13+age2225+age2629+age3035+age3644+age4554 if sex==1
replace x=hsorged+black+hispanic+married+wkless13+afdc+age2225+age2629+age3035+age3644+age4554 if sex==0

qui destring x, replace	
qui destring `cov1', replace
qui destring `cov0', replace

if `male'==1&`cov'==1{
keep if sex==1
global cov="`cov1'"
}
else if `male'==0&`cov'==1 {
keep if sex==0
global cov="`cov0'"
}
else if `male'==1&`cov'==0 {
keep if sex==1
global cov=" "
}
else {
keep if sex==0
global cov=" "
}

qui tab x, sort
qui levelsof x, local(level)
local numl0 : word count `level'
scalar numl0=`numl0'
scalar di numl0

gen flag=1
foreach l of local level {
	qui count if (x==`l' & z==0)
	scalar nullcell0=r(N)
	qui count if (x==`l' & z==1)
	scalar nullcell1=r(N)
		if (nullcell0<`mincell'|nullcell1<`mincell') {
		qui replace flag=2 if x==`l'
		}
local minl=`l'
}

qui sum flag if flag==2
	if r(N)==0 {
	*If no cells are dropped at all, need to drop one out of the total number of unique values of X
	drop if x==`minl' 
	}

qui levelsof x if flag==1, local(level)

global level= "`level'"	
local numl: word count $level
scalar numl=`numl' 
*numl is the actual number of X values used
scalar di numl
local numqt: word count $qt
scalar numqt=`numqt'
scalar dfV=numqt*numl

	set seed 123
	timer on 1
	if `bsreps'==0 {
	di in red "Need to set the number of bs replications."
	exit
	} 
	else {	
	qui test_indqt
	matrix dm=e(b)
	bootstrap _b, reps(`bsreps'): test_indqt
	matrix hatv1=e(V)
	scalar df1=dfV-diag0cnt(hatv1)
	matrix wald1=dm*invsym(hatv1)*dm'
	scalar bsstat1=wald1[1,1]
	scalar Pval1=1-chi2(df1, bsstat1)
	scalar dec1=(wald1[1,1]>invchi2(df1,0.95))
	}  	
scalar mincell=`mincell'

scalar di mincell
scalar di numl0
scalar di numl
scalar di numqt
scalar di df1
scalar di bsstat1
scalar di Pval1
scalar di dec1


timer off 1
timer list 1
scalar time=r(t1)/60
scalar di time

log close

	


